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l^j ABSTRACT 

GJ436b is a unique member of the transiting extrasolar planet population being one of the 
£> smallest and least irradiated and possessing an eccentric orbit. Because of its size, mass and 

density, GJ436b could plausibly have an atmospheric metallicity similar to Neptune (20-60 times 
solar abundances), which makes it an ideal target to study the effects of atmospheric metallicity 
on dynamics and radiative transfer in an extrasolar planetary atmosphere. We present three- 
dimensional atmospheric circulation models that include realistic non-gray radiative transfer for 
1, 3, 10, 30, and 50 times solar atmospheric metallicity cases of GJ436b. Low metallicity models 
(1 and 3 times solar) show little day/night temperature variation and strong high-latitude jets. 
In contrast, higher metallicity models (30 and 50 times solar) exhibit day /night temperature 
variations and a strong equatorial jet. Spectra and light curves produced from these simulations 
. ,-h show strong orbital phase dependencies in the 50 times solar case and negligible variations with 

^ orbital phase in the 1 times solar case. Comparisons between the predicted planet/star flux 

ratio from these models and current secondary eclipse measurements support a high metallicity 
atmosphere (30-50 times solar abundances) with disequilibrium carbon chemistry at play for 
GJ436b. Regardless of the actual atmospheric composition of GJ436b, our models serve to 
illuminate how metallicity influences the atmospheric circulation for a broad range of warm 
extrasolar planets. 

Subject headings: planets and satellites: general, planets and satellites: individual: GJ436b, methods: 
numerical, atmospheric effects 

1. Introduction 



^ETI Institute, 515 N. Wishman Road, Mountain 

View, CA 94043 The "hot Neptune" GJ436b was first discovered 



by Butler et al. (2004) and later determined to 



transit its host star as seen from earth by |Gillon 



et al. (2007b). Since the discovery of its transit 



and subsequent secondary eclipse, GJ436b has be- 



come a popular target for Hubble (Bean et al. 2008 



Pont et al.||2009[) and Spitzer flGiilon et al |: 
Deming et al.||2007| [Demory et al.||2007| |S 



2007a 



)teven- 



son et al. 2010) observations as well as modeling 
efforts ( Spiegel et al.||201Q Madhusudhan & Sea 



ger||2010[ ). Given GJ436b's mass (M p =0.0729Mj) 



and radius (i? p =0.3767Rj), its interior must con- 
tain significant quantities of heavy elements in ad- 



dition to hydrogen and helium ( Adams et al.| 2008; 
Figueira et all [2009) [Rogers &; Seager||2010[|Net-| 



telmann et al. 2010). This raises the possibil- 



ity that, like Uranus and Neptune, whose atmo- 
spheric C/H ratios lie between 20 and 60 times 
solar (iGautier et al. 1995), the atmosphere of 



GJ436b is highly enriched in heavy elements. This 
makes GJ436b an excellent case study for atmo- 
spheric chemistry, radiative transfer, and global 
circulation that should differ significantly from 
the well studied "hot Jupiters" HD209458b and 
HD189733b. 

Observations of HD 189733b using the Spitzer 
Space Telescope provided the first clear evidence 
for atmospheric circulation on an extrasolar planet 
flKnutson et al.|2QQ7||2Q09| . Most efforts to model 
atmospheric circulation for extrasolar planets have 
focused on hot Jupiters, specifically HD 189733b 
and HD209458b ([Showman fc Guillot|2002l[Show 



Cooper & S ho wman 
FTIn 2008; Menoul 



man et a l. 2008, 200 9| |Cho et al.||2003| |2008| 



2005 



20061 IDobbs-Dixon 



Rauscher 2009; Rauscher 



& Menou 2010). Only a handful of studies 



have specifically investigated the effects of non- 



synchronous rotation (Cho et al. 2008 Showman 



et al. 2009), non-zero obliquity (Langton & Laugh- 
lin 2007), and non-zero eccentricity ([Langton fc 



Laughlin 2008 ) . The possible effect of atmospheric 



composition, and hence opacity, on circulation 
patterns that may develop on extrasolar planets 
has been investigated to some extent by [Dobbs- 



Dixon & Lin (2008) and Showman et al. (2009), 



but largely ignored in most of the current two- 
dimensional and three-dimensional atmospheric 
models. Atmospheric composition is key in deter- 
mining opacity and radiative timescales that play 
a crucial role in the development of circulation on 
these planets. 

Here we present three-dimensional atmospheric 



models for GJ436b that incorporate both equi- 
librium chemistry and realistic non-gray radia- 
tive transfer. Although the actual composition of 
GJ436b's atmosphere is likely to deviate from an 
equilibrium chemistry solution as shown from the 
secondary eclipse observations of |Stevenson et alT 



(2010), our investigation still serves to explore the 



effect metallicity can play in controlling the atmo- 
spheric circulation not only on GJ436b but for a 
broad range of gaseous extrasolar planets in a sim- 
ilar temperature range. Our models are not con- 
strained to match measured chemical abundances 
and temperatures, but instead provide a system- 
atic look at how changes in atmospheric metallic- 
ity over the range from 1 times to 50 times so- 
lar values affects the basic thermal and dynami- 
cal structure of the planet's atmosphere. We do 
not expect that spectra and light curves from our 
model will provide a match to observational data, 
but instead illuminate some of the underlying at- 
mospheric physics responsible for current observa- 
tions and suggest areas of focus for future obser- 
vations. Section [2] gives an overview of the three- 
dimensional coupled radiative transfer and atmo- 
spheric dynamics model used in this study. Sec- 
tion [3] presents the global thermal structures and 
winds that develop in each of our models along 
with predicted light curves and emission spectra. 
Sections [4] and [5] provide a brief discussion of the 
results and final conclusions. 

2. Model 

The atmospheric model used in this study is a 
three-dimensional (3D) coupled radiative transfer 
and dynamics model that was specifically devel- 
oped with the study of extrasolar planetary atmo- 
spheres in mind. The Substellar and Planetary 
Atmospheric Radiation and Circulation (SPARC) 
model is described in detail in I Showman et al.l 
( [2009] ) as applied to HD189733b and HD209458b. 
A basic overview of the SPARC model along with 
the specific changes made to the model setup 
for GJ436b are presented here for completeness. 



The SPARC model employs the MITgcm (Adcroft 



et al. 2004) to treat the atmospheric dynamics 
using the primitive equations, which are valid in 
stably stratified atmospheres where the horizontal 
dimensions of the flow greatly exceed the verti- 
cal dimension. For GJ436b, the horizontal length 
scale of the flow is ~ 10 T m while the vertical scale 



height of the atmosphere is ~ 300 km. The simula- 
tions presented here take advantage of the cubed- 
sphere grid ( Adcroft et al.|2004 ) at a resolution of 
C32 (roughly 64 x 128 in latitude and longitude) to 
solve the relevant dynamic and energy equations. 
The vertical dimension in these simulations spans 
the pressure (p) range from 200 bar to 20 /ibar 
with 47 vertical levels, evenly spaced in \og(p). 
The boundary conditions in our simulations are 
an impermeable surface at the bottom and a zero 
pressure surface at the top both of which are free 
slip in horizontal velocity. 

We have coupled the MITgcm to the non- 



gray radiative transfer model of M arley fc McKay 



(1999) to realistically determine the magnitude of 
heating/cooling at each grid point. The radia- 
tive transfer model, a two-stream version of the 



Marley & McKay (1999k plane-parallel code, as- 
sumes local thermodynamic equilibrium and in- 
cludes intensities over the wavelength range from 
0.26 to 300 /im. The opacity at each pressure- 
temperature-wavelength grid point is tabulated 
using the correlated-k method (iGoody et al. 1989). 



Our extensive opacity database is described in 



Freedman et al. (2008). The chemical mixing ra- 



tios, which are computed assuming thermochem- 
ical equilibrium, are calculated as in |Lodders ~fc| 
Fegley (2002 2006). Calculated opacities assume 
a gaseous composition without particulate mat- 
ter and account for the possibility of chemical 
rainout. Because GJ436b plausibly has an atmo- 
spheric chemistry that is enhanced in heavy el- 
ements, we developed opacity tables for 3 times 
(3x), 10 times (10 x), 30 times (30 x), and 50 
times (50 x) solar metallicity in addition to the 1 
times (lx) solar metallicity opacity table. In the 
enhanced metallicity opacity tables, all elements 
other than hydrogen and helium are assumed to 
be enhanced by the same factor over current so- 
lar values. The opacity databases of Freedman 



et al.| (120081) were updated to include the opac- 
ity effects of CO2, which is an important carbon 
bearing species at higher metallicities. The full 
opacity tables are divided into 30 wavelength bins 

( [2009| . This bin- 



as outlined in Showman et al. 



ning of opacities allows for greater computational 
efficiency while only introducing small (< 1%) de- 
viations from the net radiative flux calculated with 
higher resolution opacity tables. 

For each model, the winds are assumed to ini- 



tially be zero everywhere and each column of the 
grid is assigned the same pressure-temperature 
profile. This initial pressure-temperature pro- 
file is derived from one-dimensional radiative- 
equilibrium calculations performed using the ra- 
diative transfer code in the absence of dynam- 
ics. Figure [l] shows the pressure-temperature pro- 
files derived for each metallicity case of GJ436b. 
These pressure-temperature profiles were derived 
using the methodology presented in |Fortney et al. 
(2008 2005). The physical properties assumed 



for GJ436b and its host star (GJ436A) are pre- 
sented in Table [I] Using these planetary and stel- 
lar parameters, the effective temperature (T e ff) 
of GJ436b is calculated to be 649 K, assuming 
planet-wide redistribution of the incoming stellar 
flux. This T e ff corresponds to a mean photo- 
spheric leveQ of 1 to 100 mbar depending on the 
assumed metallicity of the atmosphere (Figure pi). 
Because GJ436b is known to have an eccen- 
tric orbit, we incorporated the effects of non- 
synchronous rotation and time-varying distance 
from the host star into the SPARC model. The 
most probable rotation rate for GJ436b was de- 
termined using the following pseudo-synchronous 
rotation relationship presented in Hut] (|1981[) : 
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where P rot is the planetary rotation rate, P or 5 is 
the orbital period of the planet, and e is the eccen- 
tricity of the planetary orbit. In all cases consid- 
ered here the obliquity of the planet is assumed to 
be zero. The time-varying distance of the planet 
with respect to its host star, r(t), is determined 
using Kepler's equation ( Murray fc Dermott|1999 ) 
and used to update the incident flux on the planet 
at each radiative timestep. A diagram of GJ43b's 
orbit is presented in Figure |2j To test the impact 
of pseudo-synchronous rotation and time- varying 
stellar insolation, additional simulations for the 
1 x and 30 x solar metallicity cases were performed 
assuming synchronous rotation and zero eccentric- 
ity. 

In our models, for computational efficiency, the 
radiative timestep used to update the radiative 



1 Defined in this context as the atmospheric pressure where 
the local temperature equals the effective temperature. 



fluxes is longer than the timestep used to up- 
date the dynamics. Generally, as we increased 
the metallicity of the atmosphere, progressively 
shorter radiative and dynamical timesteps were 
needed to maintain stability. For the lx and 3x 
solar metallicity cases a dynamic timestep of 25 s 
and a radiative timestep of 200 s were used. The 
10 x and 30 x solar metallicity cases required a dy- 
namic timestep of 20 s and a radiative timestep 
of 100 s while the 50 x solar case required a dy- 
namic timestep of 15 s and a radiative timestep of 
60 s. Timestepping in our simulations is accom- 
plished through a third-order Adams-Bashforth 
scheme ( Durran||l991 ). We applied a fourth-order 
Shapiro filter in the horizontal direction to both 
velocity components and the potential tempera- 
ture over a timescale equivalent to twice the dy- 
namical timestep in order to reduce small scale 
grid noise while minimally affecting the physical 
structure of the wind and temperature fields at 
the large scale. 

We integrated each of our models until the ve- 
locities reached a stable configuration. Figure [3] 
show the root mean square (RMS) velocity as a 
function of pressure and simulated time, calcu- 
lated according to: 



Vrms(p) 



f(u 2 +v 2 )dA 
A 



(2) 



where the integral is a global (horizontal) inte- 
gral over the globe, A is the horizontal area of the 
globe, u is the east-west wind speed, and v is the 
north-south wind speed. The high-frequency vari- 
ations in the RMS velocity seen in the upper lev- 
els of both the 1 x and 50 x solar cases are largely 
due to variation in the incident stellar flux asso- 
ciated with the eccentric orbit of GJ436b. Notice 
that, in the observable atmosphere (pressures less 
than 100 mbar), the orbit- averaged winds become 
essentially steady within ~2500 Earth days for so- 
lar metallicity and ~1000 Earth days for 50 x so- 
lar metallicity. RMS wind speeds typically reach 
~1 km s _1 at photosphere levels. Any further in- 
creases in wind speeds will be small and confined 
to pressure well below the mean photosphere so 
as not to affect any synthetic observations derived 
from our simulations. As outlined in Showman 



et al. (2009) the energy available for the produc- 



tion of winds is limited largely by the global avail- 
able potential energy within the atmosphere and 



to some extent energy losses due to the Shapiro 
filter which acts as a hyperviscosity. A full discus- 
sion of the energetics of our simulated G J436b-like 
atmosphere is left for a future paper. 

3. Results 

The following sections overview the key results 
from the study of GJ436b's atmospheric circula- 
tion at lx, 3x, 10 x, 30 x, and 50 x solar metallic- 
ity. Both the thermal structure and winds in these 
simulations have a strong dependence on the as- 
sumed composition of the atmosphere for GJ436b. 
Additionally, theoretical light curves and spectra 
are produced from our 3D model atmospheres and 
compared with available data. 

3.1. Thermal Structure and Winds: De- 
pendence on Metallicity 

Figure [4] presents snapshots of the temperature 
and wind fields at three pressure levels in the at- 
mosphere for the lx and 30 x solar cases near 
secondary eclipse when the full day-side of the 
planet faces Earth (Figure [2]). Overall, the 30 x 
solar case is significantly (~ 100 K) warmer than 
the lx solar case at each pressure. The increased 
atmospheric opacity that comes with metallicity 
enhancements leads to an upward shift in the 
pressure-temperature profiles. This effect is self- 
consistently generated in the three-dimensional 
model integrations but can also be seen in the one- 
dimensional radiative-equilibrium solutions shown 
in Figure [I] Overall, the day/night temperature 
contrast in the upper layers of the atmosphere (~1 
mbar) and the equator/pole temperature contrast 
deeper in the atmosphere (~30 mbar) increase 
with atmospheric metallicity. However, because 
the pressure at a given optical depth is smaller 
at high metallicity than low metallicity, the re- 
gions that develop significant day/night tempera- 
ture contrasts shifts to higher altitude as metal- 
licity increases. At the 1 bar level, the equa- 
tor/pole temperature contrast in the 30 x solar 
case is smaller than that in the 1 x solar case. This 
lack of a strong temperature contrast at 1 bar in 
the 30 x solar case occurs because this pressure is 
at a greater optical depth in the 30 x solar case 
than in the lx solar case, and thus occurs below 
the levels with the strongest heating/cooling. 

It is also informative to compare the flow pat- 



terns indicated by the arrows in Figure [4] between 
the lx and 30 x solar cases. The overriding fea- 
ture in all simulations is the development of a pro- 
grade (eastward) flow at low pressure. The flow 
patterns in the 30 x solar case exhibit clear wave- 
like structures outside of the equatorial region at 
the 1 bar, 30 mbar, and 1 mbar levels. At the 1 
bar level, the flow is predominately westward in 
the 30 x solar case and predominately eastward in 
the lx solar case. In both the lx and 30 x solar 
metallicity cases the strength of the winds indi- 
cated by the length of the wind vectors in Figure [4] 
is a stronger function of latitude than longitude, 
except at the highest levels (1 mbar) in the 30 x 
solar case, which shows a significant day /night 
temperature contrast similar to what was seen in 
the simulations of the tidally locked hot Jupiters 
HD189733b and HD209458b from lShowman et al.l 



(2009). 



We find that the atmospheric metallicity plays 
a key role in determining the jet structure for 
a planet with temperatures similar to those ex- 
pected on GJ436b. This is demonstrated clearly 
in Figure [5j which shows the zonal-mean zonal 
wincjj versus latitude and pressure for each of the 
five atmospheric metallicities for GJ436b consid- 
ered in this study. In the lx solar case, strong 
high-latitude jets develop in the atmosphere with 
a weaker equatorial jet. Increasing the metallic- 
ity of the planet to 3 x solar causes a strengthen- 
ing of the equatorial jet and a weakening of the 
high-latitude jets. Once the metallicity of the at- 
mosphere is increased to 10 x solar or more, the 
high-latitude jets disappear and the equatorial jet 
becomes dominant. Overall, the maximum zonal 
wind speed increases with metallicity from roughly 
1300 m s _1 in the lx solar case to over 2000 m 
s _1 in the 50 x solar case. The flow is subsonic 
everywhere in our 1, 3, and 10 x solar metallic- 
ity cases. In our 30 and 50 x solar metallicity 
cases, the winds are subsonic throughout most 
of the domain, but they become marginally su- 
personic at the very top of the domain (at pres- 
sures ~1 mbar) in the equatorial region. Hydraulic 
jumps similar to those seen in the HD189733b 
and HD209458b cases presented in Showman et al.| 
( 2009 ) are present in the supersonic regions of the 



atmosphere in the 30 and 50 x solar metallicity 
cases (see 1 mbar level of the 30 x solar case in Fig- 
urepl). It is important to note that the flow in all of 
the metallicity cases are predominately eastward 
at pressure less than ~1 bar and predominately 
westward at pressures greater than ~1 bar. Mo- 
mentum conservation requires that the eastward 
momentum in the upper atmosphere of these sim- 
ulations comes from the deeper atmospheric layer, 
which requires the development of the mean west- 
ward flow at depth. The detailed mechanisms re- 
sponsible for this momentum and the jet pumping 
mechanisms themselves will be discussed in a fu- 
ture paper. 

In rapidly rotating atmospheres, atmospheric 
temperature gradients are linked to winds by dy- 
namical balances, so it is interesting to next exam- 
ine the atmospheric temperature structure in our 
simulations. Figure [6] shows the zonal- mean at- 
mospheric temperature versus latitude and pres- 
sure for the lx and 30 x solar cases. In both 
cases, the deepest isotherms (for temperatures ex- 
ceeding ~1200 K) are flat, but isotherms between 
600 and 1100 K are bowed upward, indicating a 
warm equator and cool poles. The relationship be- 
tween this structure and the winds can be under- 
stood with the thermal-wind equation, which re- 
lates the latitudinal temperature gradients to the 
zonal wind and its derivative with pressure: 



2u tan ( 



a 



2Qi 



du 
dhip 



a G(p 



2 That is, the longitudinally averaged east- west wind, where 
eastward is denned positive and westward negative. See 



Holton (2004) 



where u, </>, a, Q, p, R and T are the zonal wind 
speed, latitude, planetary radius, planetary rota- 
tion rate, pressure, specific gas constant, and tem- 
perature, respectively. The latitudinal tempera- 
ture gradient on the right-hand side is evaluated at 
constant pressure. This relationship derives from 
taking a vertical (pressure) derivative of the merid- 
ional momentum equation for a flow where the pre- 
dominant zonal-mean meridional momentum bal- 
ance is between the Coriolis, pressure- gradient, 
and curvature terms (called "gradient-wind" bal- 
ance; see Holton 2004, p. 65-68). From Equation 
O), one expects that, away from the equator, re- 
gions exhibiting vertical shear of the zonal wind 
must also exhibit latitudinal gradients of temper- 
ature. Comparing Figures [5] and [6] confirms that 
this is indeed the case: in the mid-latitudes, the 
1 x solar case exhibits the greatest vertical shear of 



the zonal wind in the pressure range of ~0.1 to 3 
bars (Figure|5|, and this is the same pressure range 
over which the greatest latitudinal temperature 
gradients occur (Figure [6]). For the 30 x solar case, 
in the mid-latitudes, the regions exhibiting sig- 
nificant wind shear are shifted upward, occurring 
from ~0.01 to less than 1 bar, and likewise this 
is the pressure range where significant latitudinal 
temperature gradients exist. The upward shift in 
the temperature gradients (and hence winds) at 
greater metallicity is the direct result of enhanced 
atmospheric opacities, which lead to shallower at- 



fect of eccentricity and rotation rate on the cir- 
culation by comparing our standard cases (Sec- 
tion 3.1) to cases with synchronous (rather than 



mospheric heating (Fortney et al. 
Dixon fc Lin||2QQ8[ ). 



2008 Dobbs- 



It is interesting to characterize the variations 
in wind speed that occur throughout the eccentric 
orbit due to the time- variable incident stellar flux. 
Figure [7] shows the RMS wind speed as a func- 
tion of pressure and simulated time with outputs 
every two hours for a ten Earth day period after 
the simulations had reached equilibrium. Over- 
all, wind speeds in these simulations of GJ436b 
vary with a frequency roughly equal to the orbital 
period at pressures above the mean photospheric 
leveF] (roughly 100 mbar in the lx solar case and 
10 mbar in the 50 x solar case). Wind speeds 
are fairly constant at pressures below the mean 
photosphere for each metallicity case of GJ436b. 
The vertical lines in Figure [7] indicate the time 
of periapse passage, which are followed by a peak 
in the RMS wind speeds. The variation in the 
wind speeds between periapse and apoapse is sev- 
eral times larger in the 50 x solar case compared 
with the lx solar case. The higher atmospheric 
metallicity cases show greater variability in wind 
speeds as function of orbital phase, which could 
affect light curves especially if hotter or more ec- 
centric systems are considered. 

3.2. Effect of eccentricity and rotation 
rate 



The models presented in Section [3TT] assume an 
eccentric orbit (hence time- variable stellar irradia- 
tion) and adopt the pseudo-synchronous rotation 
rate given in Table [I] Here, we explore the ef- 



3 In Figure|3J the high-frequency fluctuations appear to occur 
on periods of tens of Earth days, but this is an artifact that 
results from aliasing of the sampling frequency of 5 X 10 5 s 
with the orbital period. 



pseudo-synchronous) rotation rates and zero ec- 
centricity. Figure [8] presents the thermal struc- 
ture and winds at the 30 mbar level for the lx 
and 30 x solar cases assuming synchronous rota- 
tion (P ro t = P rb)- The top panels of Figure [8] as- 
sume the nominal eccentric orbit of GJ436b while 
the bottom panels assume a circular orbit with the 
nominal semimajor axis of GJ436b (Table HI). The 
assumption of synchronous rotation and/or a cir- 
cular orbit has little effect on the overall thermal 
structure and wind patterns that develop in these 
simulations. This is presumably because the ec- 
centricity of GJ436b's orbit is modest (e = 0.15) 
and changes in the average stellar flux and pseudo- 
synchronous rotation rate from a circularized and 
tidally locked orbit are small. 

Figure [9] presents the zonal-mean zonal winds 
for these same four cases (lx and 30 x solar metal- 
licity, with eccentricity of or 0.15, all using the 
synchronous rotation period). Overall, the jet 
structures differ little from the nominal cases. It 
is interesting to note that the eastward equatorial 
jet in the 30 x solar synchronous rotation cases has 
a maximum wind speed that is ~ 200 m s _1 faster 
than what is seen in the non-synchronous case. 
Because GJ436b has a relatively small eccentric- 
ity orbit the effects of non-synchronous rotation 
and time-variable heating are only small pertur- 
bations on the synchronous rotation and circular 
orbit cases. Planets with higher eccentricities are 
likely to show a larger variation in the circulation 
patterns that develop compared with circularized 
and synchronous cases. 

3.3. Light Curves and Spectra 

The SPARC model is uniquely equipped to 
produce both theoretical light curves and spec- 
tra that account not only for radiative effects, 
but also dynamic movement in the atmosphere. 
Once each of the GJ436b models reached an equi- 
librium state, pressure and temperature profiles 
were recorded along each grid column at many 
points along the planet's orbit (Figure |2j. These 
pressure-temperature profiles were then used in 
high resolution spectral calculations to determine 
the emergent flux from each point on the planet 
and which portion of that emergent flux would be 



directed toward an earth observer including limb 
darkening/brightening effects. Spectra and light 
curve generation methods are fully described in 
Fortney et al.l (|2006|). Figure fTol shows the the- 



oretical light curves, expressed as the planet/star 
flux ratio, for the 1 x and 50 x solar cases as a func- 
tion of orbital position. The 1 x solar light curves 
are relatively fiat due to the lack of a day/night 
temperature contrast as seen in Figure [4] The 
day/night temperature contrast is more prominent 
in the 50 x solar case, which results in an increase 
in the planet/star flux ratio at secondary eclipse. 
The open and filled circles in Figure [10] represent 
output from three consecutive orbits separated by 
100 simulated days to test for temporal variability 
in the light curves. Little variation is seen in the 
predicted light curves from orbit to orbit or over 
longer timescales. 

Computed emission spectra from the points 
along the orbit shown in Figure [2] are presented in 
Figure [TT] for both the lx and 50 x solar cases. As 
with all highly irradiated planets, the expectation 
is that the infrared spectra are carved predom- 
inantly by the absorption bands of H2O vapor. 
These bands are most prominent in the near in- 
frared between the J-, H-, and K-band flux peaks 
(it is at the flux peaks that the H 2 opacity is 
low). At wavelengths where H2O opacity is low, 
the deeper, hotter layers of the atmosphere can 
be seen and the emitted flux is generally higher. 
However, other molecules can imprint absorption 
features on these emission peaks, and lessen them. 
This is true in the 4-5 /im range where CO ab- 
sorbs over the redder half of this range (~4.5-5 
/im) along with CO2 (~4.3 /im), which is promi- 
nent at high metallicity. Given the assumption 
of chemical equilibrium, CH4 is abundant in all 
of the metallicity cases, which leads to a strong 
absorption band centered on 3.3 /im, as well as 
a broad band from ~7-8.5 /im. However, as the 
atmospheric metallicity is increased, the CO abun- 
dance increases linearly, and the CO2 abundance 



increases quadratically (Lodders & Fegley 2002 



Zahnle et al. 2009b), which leads to a weakening 



of the CH4 bands and a strengthening of the CO 
and CO2 bands. In the 50 x solar case a strong 
CO2 band at 4.4 /im is clearly visible that does 
not appear in the lx solar case. At high metal- 
licity, this CO2 absorption band forces flux out 
a longer and shorter wavelengths, away from the 



~4-5 /im flux peak. 

In addition to Figure 11 showing emission spec- 
tra from 1 to 30 /im, emitted flux distributions as 
a function of wavelength, with orbital phase, can 



also be analyzed. Figure [12] presents the planetary 
flux per unit wavelength for the lx and 50 x so- 
lar cases. In both the lx and 50 x solar models, 
the peak energy output of the planet occurs be- 
tween 4 and 5 /im. For the 50 x model, the strong 
CO and CO2 bands just shortward of 5 /im dra- 
matically decrease the energy output there, forc- 
ing much of the energy out at both longer and 
shorter wavelengths. The dashed line in Figure 12 
shows the integrated flux as a function of wave- 
length for secondary eclipse (blue dashed line, Fig- 
ure 12). The flux from the planet in the 1-15 /im 



range accounts for 95% of the planet's emergent 
energy, which makes this an especially important 
wavelength range for determining the atmospheric 
properties of GJ436b. 

It is interesting to note the increased variability 
in the emitted flux from the planet as a function 
of orbital position in the 50 x solar case compared 
with the 1 x solar case in Figures [TT] and [12] The 
flux emitted from the 50 x solar case at secondary 



eclipse (blue line, Figure 11) is lacking in many 
of the predominant spectral features seen at other 
orbital phases, due to a shallower day-side tem- 
perature gradient. The absorption features due to 
CH4 are much weaker at secondary eclipse indicat- 
ing a reduction in CH4 abundance on the day-side 
compared to the night-side which is seen during 
transit (orange line, Figure 11). Observing the 
flux emitted from GJ436b as a function of wave- 
length at several different points along its orbit 
could reveal a great deal about its overall chemi- 
cal composition. 

4. Discussion 

The atmospheric models presented here are 
not only useful for exploring circulation regimes, 
chemistry, and radiative transfer, but can also 
provide insight into current observations and 
help guide future observations of GJ436b. Fig- 
ure[l0]also includes the available Spitzer secondary 
eclipse measurements from | Stevenson et" al. (2010). 
In all metallicity cases, our predicted planet/star 
flux ratio falls short of the measured 3.6, 5.8, 8.0, 
16.0, and 24.0 /im values and is higher than the 



observed 4.5 /im value. However, it is useful to 
note that 50 x solar model planet/star flux ratio 
comes much closer to matching the observations 
than the lx solar model, which could hint at a 



high metallicity as already suggested by|Stevenson 



et al.| ([2010]). High metall icity solutions ( ~30 x 
solar) are also favored by Spiege fet al.| ( |2010 ) 
to match the 8 /mi observations of GJ436b from 



Deming et al. (2007). The predicted planet/star 



flux ratios in the 50 x solar case are within two 
sigma of the values measured by | Stevenson et al. 



(2010) at all bandpasses except 3.6 /im. 

In comparing our predicted planet/star flux ra- 
tios with those observed by Stevenson et aT] (2010) 
it is important to remember that we have not al- 
tered the temperature or chemistry in our mod- 
els in an attempt to match observations. The in- 
terplay between the equilibrium chemistry mixing 
ratios, the absorption and emission of flux, and 
the atmospheric dynamics dictates the tempera- 
ture structure and emergent spectrum. As pointed 



out in Stevenson et al. (2010) it is likely that dis- 



equilibrium chemistry plays a strong role in the 
atmosphere of GJ436b. They suggest a CO/CH4 
ratio that is many orders of magnitude larger than 
what one would predict from an equilibrium chem- 
istry model. Enhancement of CO at the expense 
of CH4 is well known in the atmospheres of the so- 
lar system's giant planets and brown dwarfs, but 
not to such an extreme degree. Additionally, they 
suggest that the photochemical destruction of CH4 
is important to further lessen this molecule's im- 
portance in the planet's atmosphere. The carbon 
chemistry of an atmosphere will strongly affect 
flux measurements in the 3.6, 4.5, 5.8, and 8.0 



/mi bandpasses. As shown in Figure 11 and dis- 
cussed in Section [331 lowering the amount of CH4 
in the atmosphere will increase the emergent flux 
from the planet in the 3.6 and 8.0 /im bandpasses. 
Higher order hydrocarbons produced by mixing 
and photochemistry ( Zahnle et al.|2009a ), such as 
C2H2, C2H4, and C2H6, are also strong absorbers 
throughout the near- and mid-infrared. Increasing 
the amount of CO2, along with CO, in the atmo- 
sphere will decrease the emergent flux from the 
planet in the 4.5 /im bandpass while at the same 
time increasing the emergent flux at wavelength 
on either side of this bandpass, including the 3.6 
and 5.8 /im bandpasses. 

Although our models do not include disequi- 



librium chemistry, they can provide an important 
constraint for disequilibrium chemistry models — 
namely, estimates of dynamical timescales and 
vertical mixing rates. Disequilibrium chemistry is 
expected in all regions of the atmosphere where 
dynamic timescales, Tdyn-, are shorter than chem- 
ical timescales, Tchem- The dynamic timescale is 
given simply by 



7~dyr t 



L 

v 



(4) 



where L is the relevant length scale in the hor- 
izontal or vertical direction and V is horizontal 
or vertical wind speed. For one-dimensional pho- 
tochemical models, the timescales considered are 
only those in the vertical direction in which case 
V — w, where w is the vertical velocity as a 
function of height or pressure in the atmosphere. 
Given the values of uo = Dp/ ' Dt from our models, 
the vertical velocity can be estimated as 



-H 

c 

P 



(5) 



where H is the scale height of the atmosphere, 
which is a function of p the pressure at given at- 
mospheric level. Figure [13] shows the RMS values 
for w calculated as global, day-side, and night-side 
averages for both the lx and 50 x solar cases at 
secondary eclipse. The RMS values for w were cal- 
culated as wrms(p) — \l A~ x 1 w<2 dA, where the 
integral is a horizontal integral (at constant pres- 
sure) over the day-side, night-side, or entire globe 
as appropriate. In both the lx and 50 x solar 
cases the vertical wind speeds increase monotoni- 
cally from pressures around 1 bar to 0.1 mbar with 
peak speeds around 22ms" 1 . The vertical eddy 
diffusion coefficient, K zz , can be calculated from 
these vertical wind speed profiles simply by mul- 
tiplying by the relevant vertical length scale, L. 
As shown in iSmithl ([1998 ) , in order to calculate 



the correct value of L at each pressure level both 
dynamic and chemical timescales must be consid- 
ered, but to first order L = H especially near the 
quench level where Td yn — ^chem- Since H ~ 300 
km for GJ436b, one can expect K zz values from 
10 8 cm 2 s _1 near 100 bar to 10 11 cm 2 s _1 near 0.1 
mbar. These values of K zz are in line with those 



favored by Madhusudhan & Seager (2010) to ex- 



plain the disequilibrium CO/CH4 ratio needed to 
fit the Stevenson et al. (|2010) observations. 



Regardless of the detailed chemistry, our mod- 
els suggest that the basic circulation regime of 
planets in the temperature range of GJ436b de- 
pend strongly on overall metallicity. High metal- 
licity cases (30-50 x solar) produce a dominant 
eastward equatorial jet that peaks on or near the 
equator (Figure [5]). This pattern resembles that 
previously obtained for more strongly irradiated 
hot J up iters (Cooper fc Showman|2005^ 



et al. 2008, 2009 Rauscher & Menou 



Showman 



2010). At 



low metallicity, on the other hand, the equato- 
rial jet weakens and fast eastward jets develop 
in midlatitudes — a jet pattern distinct from those 
previously reported in the hot- Jupiter modeling 
literature. These changes suggest that the atmo- 
sphere experiences a regime shift where different 
dynamical mechanisms control the jet structure at 
low versus high metallicity. It may be that each 
regime is relevant to a range of planets of differing 
effective temperatures and compositions, provid- 
ing a motivation to better understand their un- 
derlying dynamics. We will discuss the specific 
mechanisms responsible for this regime shift in a 
future paper. 

The atmospheric circulation models presented 
here represent a first step in better understand- 
ing the dynamical, radiative, and chemical mech- 
anisms that shape the atmosphere of GJ436b. Al- 
though chemical disequilibrium is likely to play 
an important role, it is important to first consider 
— as we have done here — a baseline atmospheric 
model where the effects of chemical equilibrium 
opacity are tested before more complex chemistry, 
clouds, and other factors are added and compli- 
cate the interpretation of the results. It is unlikely 
that the presence of disequilibrium chemistry in 
the atmosphere of GJ436b would strongly affect 
the global circulation patterns seen in this study. 
|Madhusudhan fc Seager| ( |2010[ ) suggest that a 
high metallicity (~30x solar) atmosphere, with 
chemical disequilibrium induced both by thermal 
quenching and photochemistry, can best explain 



the observations of Stevenson et al. (2010). If this 



is the case, then one would expect global circula- 
tion patterns on GJ436b akin to our 30 x and 50 x 
solar models with a strong equatorial jet and non- 
negligible day /night temperature contrasts. As 
shown in Figures [10] and [TT] one would expect to see 
variations in the planet/star flux ratio as a func- 
tion of orbital phase for a high metallicity G J436b- 



like atmosphere. This variation in the planet/star 
flux ratio could reveal a great deal about circu- 
lation patterns on GJ436b along with day/night 
chemistry differences. For this reason we recom- 
mend that full-orbit light curves of the planet in 
the 3.6 /im bandpass be obtained during the Warm 
Spitzer mission. Looking further into the future, 
with the James Webb Space Telescope (JWST) a 
combination of NIRSpec or NIRCam data from 3- 
5 /im, as well as MIRI data shortward of ~15 /im, 
will sample the vast majority of the planet's emit- 
ted energy. Observations as a function of orbital 
phase with these instruments could tightly con- 
strain the energy budget of the planet. Detection 
of the CO2 band depth at 4.4 /im with NIRSpec 
would be also an important metallicity indicator 
and probe of the carbon chemistry of GJ436b. 

5. Conclusions 

In this study we have investigated atmospheric 
circulation for the hot Neptune GJ436b for various 
assumed atmospheric metallicities using a three- 
dimensional coupled radiative transfer and general 
circulation model. We have found that assumed 
atmospheric composition for a planet like GJ436b 
can have strong effects on both day/night recir- 
culation and jet patterns within the atmosphere. 
Light curves and spectra produced from our mod- 
els show that enhancements in atmospheric metal- 
licity over solar produce an increase as well as or- 
bital phase variations in the planet/star flux ratio. 
Given the possible strong dependence of emitted 
flux with orbital phase, GJ436b could provide an 
important probe of atmospheric chemistry outside 
of our solar system. Moreover, we showed that for 
warm gaseous extrasolar planets in the effective 
temperature range of ^600-900 K, there exists a 
regime shift as metallicity is increased from a cir- 
culation dominated by mid-latitude jets and min- 
imal longitudinal temperature differences to one 
dominated by an equatorial jet and large day- night 
temperature differences. Although these models 
of GJ436b do not include all the processes at play 
in any given atmosphere, the basic trends for at- 
mospheric circulation as a function of metallicity 
will be useful in better understanding GJ436b and 
many other extrasolar planetary atmospheres. 
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Fig. 1. — Initial pressure temperature profiles assumed for each metallicity case of GJ436b. Lines of equal 
abundance for CH4 vs CO and NH3 vs N2 are shown to highlight the dominant carbon and nitrogen bearing 
species at each pressure level for each metallicity case. As the metallicity in the atmosphere is increased form 
lx to 50 x solar, the dominant carbon bearing species changes from CH4 to CO. Diamonds represent the 
level of the mean photosphere (T = T e ff), which decreases in pressure as the metallicity is increased. Profiles 
assume planet-wide redistribution of absorbed incident energy, (a color version of this figure is available in 
the online journal.) 
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Fig. 2. — Orbit of GJ436b. The true anomaly, /, represents the angular distance of the planet from periapse. 
Assuming a longitude of pericenter, w, of 343° from |Deming et aL (2007), transit occurs at / = 107° and 



secondary eclipse occurs at / = —73°. Dots along the orbital path represent points where data was extracted 
to produce Figures 10, 11, and 12 Colored dots represent points near periapse (red), transit (orange), 
apoapse (green), and secondary eclipse (blue), which correspond to the colored spectra presented in Figures 
11 and [12] Figure is to scale with the small purple dot after periapse representing the size of G J436b in 
relation to its host star and orbit, (a color version of this figure is available in the online journal.) 
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Fig. 3. — RMS velocity (colorscale) as a function of pressure and simulated time for the lx (top) and 50 x 
(bottom) solar cases of GJ436b. The RMS velocity at each pressure level is calculated from the instantaneous 
wind speeds recorded every 5xl0 5 s. Simulations are continued until the winds reach a relatively flat profile 
at all pressure levels, (a color version of this figure is available in the online journal.) 
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Fig. 4. — Temperature (colorscale) and winds (arrows) for the lx (left) and 30 x (right) solar metallicity 
cases of GJ436b. For both lx and 30 x solar cases the thermal structure and winds are shown at the 1 
mbar (top), 30 mbar (middle), and 1 bar (bottom) levels of the simulation. The longitude of the substellar 
point is indicated by the solid vertical line in each panel. Each panel is a snap shot of the atmospheric state 
taken near secondary eclipse (/ = —73°, Figure [2]). The horizontal resolution of these runs is C32 (roughly 
128x64 in longitude and latitude) with 47 vertical layers, (a color version of this figure is available in the 
online journal.) 
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Fig. 5. — Zonal- mean zonal winds for the five atmospheric metallicities considered in this study for GJ436b 
assuming pseudo-synchronous rotation. The wind speeds presented here represent 100 day averages of the 
zonal winds taken after each simulation was considered to have reached an equilibrium state. The colorbar 
shows the strength of the zonally averaged winds in m s _1 . Contours are spaced by 100 m s _1 . Positive 
wind speeds are eastward, while negative wind speeds are westward. Note the significant change in the jet 
structure as a function of atmospheric metallicity. (a color version of this figure is available in the online 
journal.) 
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Fig. 6. — Zonal-mean temperatures (colorscale) as a function of pressure and latitude for the pseudo- 
synchronous lx (top) and 30 x (bottom) solar cases of GJ436b. Contours represent isotherms and are 
spaced by 50 K. Gradients in temperature along an isobaric surface tend to drive zonal winds according 
to equation |3|. The 30 x solar case shows stronger thermal gradients at lower pressures that extend into 
lower latitudes when compared with the lx solar case. This change in the thermal gradient structure can 
be directly related to the zonal wind profiles seen in Figure [5j (a color version of this figure is available in the 
online journal.) 
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Fig. 7. — RMS velocity (colorscale) as a function of pressure and simulated time for the 1 x (top) and 50 x 
(bottom) solar cases of GJ436b as calculated from high- cadence simulation outputs every 7200 s. Vertical 
lines indicate periapse passage. Overall, wind speeds in the atmosphere vary with a period equal to the 
orbital period. The peak in the wind speeds typically occurs 4-8 hours after periapse passage, (a color 
version of this figure is available in the online journal.) 
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Fig. 8. — Temperature (colorscale) and winds (arrows) for the lx (left) and 30 x (right) solar metallicity 
cases of GJ436b assuming synchronous rotation at the 30 mbar level. The top panels represent simulations 
where synchronous rotation was assumed, but the nominal eccentric orbit was maintained. The bottom 
panels represent simulations where synchronous rotation and a circular orbit with the nominal semimajor 
axis (Table [I]) were assumed. In the eccentric cases (top) the panels represent a snap shot taken near 
secondary eclipse (/ = —73°, Figure [2]). The solid vertical line in each panel represents the longitude of the 
substellar point. (a color version of this figure is available in the online journal.) 
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Fig. 9. — Zonal-mean zonal winds for the synchronous rotation (P ro t — Porb) cases at lx (left) and 30 x 
(right) solar metallicity GJ436b atmospheres. The top panels represent the jet structure for a synchronously 
rotating GJ436b in an eccentric orbit, while the bottom panels assume a circular orbit with the same 
semimajor axis. The wind speeds presented here represent 100 day averages of the zonal winds taken 
after each simulation was considered to have reached an equilibrium state. Note that jet structures for the 
synchronous cases are very similar to the zonal-mean zonal wind plots presented in Figure [5] for the same 
atmospheric metallicity. (a color version of this figure is available in the online journal.) 
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Fig. 10. — Planet/Star flux ratio as a function of orbital phase in each of the Spitzer bandpasses for the lx 
and 50 x solar metallicity cases of GJ436b. The mean anomaly is an angle that increases linearly with time 

o 

from periapse passage. For GJ436b, transit and secondary eclipse occur at a mean anomalies of 90 and 

o 

303 respectively. The filled and open circles represent data extracted with 100 simulated days of separation 
to test for temporal variability, which appears to be minimal. The secondary eclipse measurements from 



Stevenson et al. (2010) are shown for comparison, (a color version of this figure is available in the online 
journal.) 
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Fig. 11. — Flux per unit frequency, F v (erg s _1 cm -2 Hz -1 ), as a function of wavelength for the lx (top) 
and 50 x (bottom) solar metallicity cases of GJ436b. The black spectra presented for each case are taken 
from 32 locations along a single orbit as shown in Figure [2] The central wavelengths of J-, H-, K-, L-, and 
M-bandpasses are indicted by the corresponding letters at the top of the plot. Dotted lines at the bottom 
indicate the bandpasses of the four Spitzer IRAC bands from 3-9 /im, IRS blue filter at 16 /im, and the MIPS 
24 /im bandpass. Colored spectra are taken from secondary eclipse (blue), periapse (red), transit (orange), 
and apoapse (green). The 50 x solar case shows a strong dependence of the emergent flux density on orbital 
position while the lx solar case does not. (a color version of this figure is available in the online journal.) 
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Fig. 12. — Flux per unit wavelength, F\ (erg s _1 cm -2 /im _1 ), as a function of wavelength for the lx (top) 
and 50 x (bottom) solar metallicity cases of GJ436b. The spectra presented for each case are taken from 
several points along the orbit as shown in Figure [2j with secondary eclipse (blue), periapse (red), transit 
(orange), and apoapse (green) highlighted. The dashed line represents the integrated flux as a function of 
wavelength from the planet at secondary eclipse and is tied to the axes on the right with 1.0 indicating the 
total integrated flux over all wavelengths. Dotted lines at the bottom indicate the bandpasses of the four 
S^itzer IRAC bands from 3-9 /im. Solid gray bars at the top indicate the wavelength coverage of 3 planned 
instruments for the JWST: NIRCam, NIRSpec, and MIRI. The vertical position of the bars is arbitrary and 
does not signify instrument sensitivity. (a color version of this figure is available in the online journal.) 
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Fig. 13. — RMS vertical velocity (wrms) as a function of pressure for both the lx (left) and 50 x solar 
metallicity cases of GJ436b at secondary eclipse. Global (green line), day-side (red line), and night-side 
(blue line) averages of the vertical velocity are presented. These wrms profiles are useful in determining 
the vertical mixing rate in the atmosphere as a function of pressure for use in disequilibrium chemistry and 
photochemical models. Note that the decrease in wrms near the top of the domain results from the existence 
of the model's upper boundary where w is forced to be zero. The vertical velocities would likely continue 
increasing with altitude if the top of the model had been placed at even lower pressures, (a color version of 
this figure is available in the online journal.) 
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Table 1: GJ436A/b parameters. 



Parameter Value 



R p (Rj) 


0.3767 


M p (Mj) 


0.0729 


g (m s~ 2 ) 


12.79 


a (AU) 


0.02872 


e 


0.15 


zu (deg) 


343 


P or6 (days) 


2.64385 


P rot (days) 


2.32851 


i?* (R ) 


0.464 


M, (M ) 


0.452 


Teff (K) 


3350 



Note. — Planetary and stellar parameters taken from plbrres et aT7| ( |2008| >. Values for e and vo were taken from |Deming et al.| 
p007] ). 
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